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We present results of recent high-statistics Monte Carlo simulations of the Edwards-Anderson 
Ising spin-glass model in thiee and four dimensions. The study is based on a non-Boltzmann 
sampling technique, the multi-overlap algorithm which is specifically tailored for sampling 
rare-event states. We thus concentrate on those properties which are difficult to obtain with 
standard canonical Boltzmann sampling such as the free-energy barriers Fg in the probability 
density Pj{q) of the Parisi overlap pai'ameter q and the behaviour of the tails of the disorder 
averaged density P{q) = [Pj{q)]a.v- 



1 Introduction 

A widely studied class of spin-glass materials^^^"'^^^ consists of dilute solutions of mag- 
netic transition metal impurities in noble metal hosts, for instance^ Au-2.98% Mn or^ 
Cu-0.9% Mn, which is one of the best investigated metallic spin glasses. In these systems, 
the interaction between impurity moments is caused by the polarization of the surround- 
ing Fermi sea of the host conduction electrons, leading to an effective interaction of the 
so-called RKKY form^ 

cosM^ (1) 

where kp is the Fermi wave number. For an illustration, see Fig. 1. This constitutes the 
two basic ingredients necessary for spin-glass behaviour, namely 



randomness - in course of the dilution process the positions of the impurity moments 
are randomly distributed, and 

competing interactions - due to the oscillations in (1) as a function of the distance R 
between the spins some of the interactions are positive and some are negative. 



The competition among the different interactions between the moments means that no 
single configuration of spins is uniquely favoured by all of the interactions, a phenomenon 
which is commonly called "frustration". This leads to a rugged free-energy landscape 
with probable regions (low free energy) separated by rare-event states (high free energy), 
illustrated in many previous articles by vague sketches similar to our Fig. 2. Experimentally 
this may be inferred from the phenomenon of aging which is typical of measurements of 
the remanent magnetization in the spin-glass phase^. 
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Figure 1. The two basic ingredients for spin-glass behaviour: Randomly distributed spin moments of transition 
metal impurities (e.g. Mn) in a noble metal host (e.g. Cu), and the characteristic oscillatory form of their effective 
RKKY interactions with competing positive (ferromagnetic) and negative (antifen'omagnetic) regions. 

However, despite the large amount of experimental, theoretical and simulational work 
done in the past thirty years to elucidate the spin-glass phase, ^'^ •^ '^ the physical mech- 
anisms underlying its peculiar properties are not yet fully understood. To cope with the 
enormous complexity of the problem various levels of simplified models have been stud- 
ied theoretically. A simplified lattice model which reflects the two basic ingredients for 
spin-glass behaviour is the Edwards-Anderson^ Ising (EAI) model defined through the 
energy 

E ^ - ^Jik SiSk , (2) 

where the fluctuating spins Si can take the values ±1. The coupling constants Jik are 
quenched, random variables taking positive and negative signs, thereby leading to com- 
peting interactions. In our study we worked with a bimodal distribution, Jn, — ±1 
with equal probabilities, but other choices such as the Gaussian distribution, V{Jik) oc 
exp(— Jj^^/2A^) with A parameterizing its width, have also been considered, in particu- 
lar in analytical work. In (2), the lattice sum runs over all nearest-neighbour pairs of a 
d-dimensional (hyper-) cubic lattice of size N = with periodic boundary conditions. 

An analytically more tractable mean-field model, commonly known as the Sherrington- 
Kirkpatrick^" (SK) model, emerges when each spin is allowed to interact with all others. 
Alternatively one may consider the mean-field treatment as an approximation which is 
expected to become accurate in high dimensions'^. In physical dimensions, however, its 
status is still unclear and an alternative droplet approximation'^ has been proposed. The 
two treatments yield conflicting predictions which have prompted quite a controversial 
discussion over many years. Numerical approaches such as Monte Carlo (MC) simulations 
can, in principle, provide arbitrarily precise results in physical dimensions. In practice, 
however, the simulational approach is severely hampered by an extremely slow (pseudo-) 
dynamics of the stochastic process. 
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Figure 2. Typical sketch of the rugged free-energy landscape of spin glasses with many minima separated by 
rare-event barriers. 



To overcome the slowing-down problem various novel simulation techniques have been 
devised in the past few years. While some of them only aim at improving the (pseudo-) 
dynamics of the MC process, others are in addition also well suited for a quantitative char- 
acterization of the free-energy barriers responsible for the slowing-down problem. Among 
the latter category is the multi-overlap algorithm^'^ which has recently been employed in 
quite extensive MC simulations^*^^^'^^^^^ of the EAI spin-glass model in three and four di- 
mensions. The purpose of this note is to give an overview of the recent progress achieved 
with this method. 



2 Model Parameters and Simulation Method 

As order parameter of the EAI model one usually takes the Parisi overlap parameter^ ^ 

i=l 

where the spin superscripts label two independent (real) replicas for the same realization 
of randomly chosen exchange coupling constants J' ~ {Jik}- For given J7 the probability 
density of q is denoted by Pj{q), and thermodynamic expectation values are computed 
as{...)j= E{.}(-..)exp(-/3i/[J])/EMexp(-/3i/[J]), where /? = l/fesT is the 
inverse temperature in natural units. The freezing temperature is known to be at /3c = 
0.90(3) in 3D (Ref. 18) and at /3c = 0.485(5) in 4D (Ref. 19), respectively. 

On finite lattices the results necessarily depend on the randomly chosen quenched dis- 
order To get a better approximation of the infinite system (apart from special problems 
with so-called non-self-averaging), one performs averages over many hundreds or even 
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thousands of (quenched) disorder realizations denoted by 



where # ^7 {-^ oo) is the number of reahzations considered. Below the freezing tem- 
perature, in the infinite- volume limit N ^ oo, a non- vanishing part of P{q) between its 
two delta-function peaks at ±qmax characterizes the mean-field picture^^ of spin glasses, 
whereas in ferromagnets as well as in the droplet picture^^ of spin glasses P{q) exhibits 
only the two delta-function peaks. Most studies so far considered mainly the averaged 
quantities. 

For a better understanding of the free-energy barriers sketched in Fig. 2, the probabil- 
ity densities for individual realizations J play the central role. It is, of course, impossible 
to get complete control over the full state space, and to give a well-defined meaning to 
"system state" (the x-axis in Fig. 2), one has to concentrate on one or a few characteris- 
tic properties. In our work we focused on the order parameter and thus those free-energy 
barriers F"^ which are reflected by the minima of Pj{q). A few typical shapes of Pj{q) 
as obtained in our simulations are shown in Fig. 3. Conventional, canonical MC simula- 
tions are not suited for such a study because the likelihood to generate the corresponding 
rare-event configurations in the Gibbs canonical ensemble is very small. This problem is 
overcome by non-Boltzmann sampling^'''^^ with the multi-overlap weight^^ 



where the two replicas are coupled by Sj{q) in such a way that a broad multi-overlap 
histogram Pj"'^{q) over the entire accessible range —1 < q < 1 is obtained, see Fig. 4. 
When simulating with the multi-overlap weight (5), canonical expectation values of any 
quantity O can be reconstructed by reweighting, (O) = {Oe^^^)j/{e~^^) j. 
The multi-overlap algorithm may be summarized as follows: 

• An iterative construction of the weight function Wj{q) = exp{Sj{q)), for each of 
the quenched disorder realizations separately. 

• An equilibration period with fixed weight function. 

• The production run with fixed weight function. 

Ideally Wj should satisfy Pj^'^{q) = Pj^'^{q)Wj{q) = const., i.e., should give 
rise to a completely flat multi-overlap probability density Pj^^'^{q) as sketched in Fig. 4. 
Of course, P^^{q) is a priori unknown and one has to proceed by iteration. Let us thus 
assume some approximate Wj_n is given. The simulation would then yield Pj^.^ which, in 
general, is not yet perfectly flat. If Pj"^^ was sampled with arbitrary precision, the desired 
weight function would be Wj oc Wj^n/Pj^n- update procedure we are actually 

only interested in relative transition amplitudes and it is therefore useful to rephrase the 
iteration in terms of 



Rj{q) ^Wj{q + ^q)/Wj{q) = RjA<i)Pj7{i)IPj7{i + - (6) 



J 



J 




(5) 
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where Aq = 2/N is the step-width in q as defined in Eq. (3). When updating the spins 
in the nth iteration, Rj n{q) — Wj,n{q + ^Q)/Wj_n{q) has to be considered as a given, 
fixed function of q. The multi-overlap histogram P'^^^{q), however, is always a noisy 

estimator whose statistical errors can be estimated by ^JPj^^isi) (neglecting auto- and 
cross-correlations and assuming unnormalized histograms counting the number of hits). 
By taking the logarithm in Eq. (6) it is then straightforward to obtain the squared error on 
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Priq)exp[Sj{q)] 







Figure 4. Illustration of the relation between canonical densities P^^" (q) depicted in Fig. 3 and ideally flat 
multi-overlap densities Pj^'^{q). 



(7) 



We now have two noisy estimators, Rj{q) and Rj,n{q) (with squared inverse er- 
rors p{q) and which may be Hnearly combined to yield an optimized estimator 

\iiRj^n+i{q) = K„{q)\nRj{q) + [1 - Kn{q)\ \nRj.a{q), with 



(q) ^ p{q)/[p{q) +Pn{q)] 



(8) 



such that the statistical error of the linear combination is minimized. By exponentiating 
the optimized estimator and using Eq. (6), we finally arrive at the recursion 



Rj,n+l(q) = RAqr-'-'^RjAqy^''''^'^ 

= RjAq) [P7„^(g)/P^:r(9 + Ag) 
Pn+i{q) =p{q) +Pn{q) ■ 



(9) 
(10) 



Once Wj{q) — exp{Sj{q)) is determined and kept fixed, the system is equilibrated and 
the data production can be performed. 

We measure the (pseudo-) dynamics of the multi-overlap algorithm in terms of the au- 
tocorrelation time Tj^'^ which is defined by counting the average number of sweeps it takes 
to complete the cycle q = ^ \q\ = 1 ^ Q. Adopting the usual terminology^°'^^ for a 
first-order phase transition, we shall call such a cycle a "tunnelling" event. The weight iter- 
ation was stopped after at least 10 "tunnelling" events occurred, and in the production runs 
we collected at least 20 "tunnelling" events. To allow for standard reweighting in tempera- 
ture we stored besides Pj{q) and the time series of q also the energies and magnetizations 
of the two replicas. The number of sweeps between measurements was adjusted by an 
adaptive data compression routine to ensure that each time series consists of 2^^ = 65 536 
measurements separated by approximately r"? ""^ sweeps. 
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3D 


4D 


L 


T = 1.00 


T = 1.14 


T = 1/0.6 


4 


8192 


8192 


4096 


6 


8192 


8192 


4096 


8 


8192 


8192 


1024 


12 


640 


1024 




16 




256 





Table 1. Number of simulated realizations #J. 



3 Results 

Due to the large number of realizations to be simulated, the final results are relatively 
costly. The studied cases are summarized in Table 1, where also the simulation tempera- 
tures are given: T = 1 « 0.88Tc and T = 1.14 « in 3D, and T = 1/0.6 w 0.85Tc in 
4D. The Jik reaUzations were drawn using the pseudo random number generators RAN- 
MAR and RANLUX ^^'^^ (luxury level 4). For the spin updates we always employed 
the faster RANMAR generator. 

By fitting the averaged autocorrelation times to the power-law ansatz ln([T™"'']av) = 
a + z In(iV), we obtainedi'"^ z = 2.32(7) and z = 1.94(2) in the 3D and 4D spin-glass 
phase, respectively. Even though the quality of the fits is quite poor and an exponential be- 
haviour can hardly be excluded, they clearly show that the slowing down is quite off from 
the theoretical optimum z = 1 one would expect if the multi-overlap autocorrelation time 
Tj'^^ is dominated by a random- walk behaviour between q = —1 and +1. In multicanon- 
ical simulations of the 3D model with broad energy histograms an even larger exponent of 
z = 2.8(1) has been observed^^. The large values of z suggest that the canonical overlap 
barriers are not the exclusive cause for the slowing down of spin-glass dynamics below the 
freezing point, i.e., the projection of the multi-dimensional state space onto the g-direction 
hides important features of the free-energy landscape of the model. 



3.1 Free-Energy Barriers 



To define effective free-energy barriers we first construct an auxiliary ID MetropoUs- 
Markov chain^^ which has the canonical Pj{q) probability density as its equilibrium dis- 
tribution. The tridiagonal transition matrix 



/I - ^2,1 
W2,l 1 



T = 



Wl,2 







W2,3 
W2,3 - W4,3 
"^4,3 



V 



(11) 



is given in terms of the probabilities Wij = \ min^l, Pjiqi)/Pj{qj)j {i 7^ j) for jumps 

from state q = Qj to q = qt = i/N in steps of Aq = ±2/N or 0. Since T fulfills the 
detailed balance condition (with Pj) it has only real eigenvalues. The largest eigenvalue 
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Figure 5. Distribution function Fq (14) for the 3D overlap barriers (13) (left) and the energy (right) in the spin- 
glass phase at T = 1 in units of their median values. 



Aq equals unity and is non-degenerate. The second largest eigenvalue Ai determines the 
autocorrelation time of the chain (in units of sweeps), 



TVlnAi 7V(1-Ai) ' 



(12) 



which we use to define the associated effective free-energy barrier in the overlap parameter 
q as 

F's = Hrl) . (13) 

Our finite-size scaling (FSS) analyses of the thus defined overlap barriers are based on 

the (cumulative) distribution function F{x). More precisely, we constructed^^ a peaked 
distribution function Fq{x) by reflecting F{x) at its median value 0.5, 

^Qy^> -\l-F{x) iorF{x)>0.5. ^ ^ 

For self-averaging data the function Fq collapses in the infinite- volume limit to Fq (x) = 
0.5 for X — [^]av 3rid otherwise. For non-self-averaging quantities the width of Fq 
stays finite. The concept carries over to quantities which diverge in the infinite-volume 
limit, when for each lattice size scaled variables x/xmed are used, where Xmed denotes the 
median defined through F(a;mod) = FQ{xme<i) = 0.5. 

The behaviour of FQ{Fg/ Fg J shown on the l.h.s. of Fig. 5 for the 3D case at 
T — 1 clearly suggests that Fg is a non-self-averaging quantity. In 4D the evidence is even 
stronger than in 3D. Non-self-averaging was also observed^^ for the autocorrelation times 
Tj^'^ of our algorithm while the energy is an example for a self-averaging quantity; cf. the 
r.h.s. of Fig. 5. For non-self-averaging quantities one has to investigate many samples and 
should report the FSS behaviour for fixed values of the cumulative distribution function F. 

In Ref. 15 we performed FSS fits for F = i/16, i = 1, . . . , 15, assuming an ansatz 
suggested by mean-field theory^®'^^, 

F^ = ai+ a2 N^'^ , (15) 
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Figure 6. FSS fits of the 3D overlap barriers Fg in the spin-glass phase at T = 1 for fixed values of the 
distribution function, F = i/16, i = 1, . . . , 15 (from bottom to top). Shown are the results for the mean-field 
prediction (15) (left) and the logarithmic ansatz (16) (right). 



corresponding to oc e°'^'^^'^ . In both dimensions the goodness-of-fit parameter^° Q 
turned out to be unacceptably small. The 3D fits are depicted on the l.h.s. of Fig. 6. We 
therefore also tried fits to the ansatz 

Fl = ln(c) + a ln{N) , (16) 

corresponding to oc N"; cf. the r.h.s. of Fig. 6. Since in 3D as well as in 4D the average 
(5-value is now within the statistical expectation, the latter ansatz (16) is strongly favoured 
over the mean-field prediction (15). As a function of F {= 1/16 - 15/16) the exponent 
a = a{F) in the power law (16) varies smoothly from 0.8 to 1.1 in 3D and from 0.8 
to 1.3 in 4D. A similar analysis^^ for the autocorrelation times r™"'' of the multi-overlap 
algorithm gives exponents a{F) which are larger, a'^^'^{F) « Q:^(F) -|- 1. This is in 
agreement with our observation that other relevant barriers exist, which cannot be detected 
in the overlap parameter q. 



3.2 Averaged Probability Densities P{q) 

The averaged canonical densities P{q) of the 3D model are shown in Fig. 7 for both T = 
1 w 0.88Tc and T = 1.14 w T^- At least close to one expects that, up to finite-size 
corrections, the probability densities scale with system size. A method to confirm this 
visually is to plot P'{q) = crP{q) versus q' = q/a, where a is the standard deviation. By 
fitting the standard deviation to the expected FSS form a = cxL-l^/" we obtained^'^, 

0.312(4), Q = 0.32 for T= 1.14, and (17) 
0.230(4) , Q = 0.99 for T = 1 . (18) 

On the l.h.s. of Fig. 8 we show the scaling plot^^ for T = 1.14 which demonstrates that 
the five probability densities collapse onto a single master curve. 
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Figure 7. Overlap probability densities for the 3D EAI model in the spin-glass phase (left) and close to criticality 
(right). 



3.3 TaUs of P{q) 

The multi-overlap algorithm becomes particularly powerful when studying the tails of the 

probability densities which are highly suppressed compared to the peak values; see the 
r.h.s. of Fig. 8 which shows P{q) at T = 1.14 over more than 150 orders of magni- 
tude. Based on the repUca mean-field approach, theoretical predictions for the scaling 

behaviour of the tails have been derived by Parisi and collaborators'^'^. They showed that 
P{q) = Pmax .f{N [q — (?j^ax)'^) 1 > 9inax concludcd more quantitatively that 

P{q) ~ cxp [-C1 N{q- for N {q - q^^^f large , (19) 

with a mean-field exponent of a; = 3. By allowing for an overall normalization factor Cq^^ 
and taking the logarithm twice we have performed fits of the form^^'^^ 

F = In [- ln(P/4^) )] - In = In ci + X ln(g - q^^^) . (20) 

Leaving the exponent a; as a free parameter, we arrived at the estimates x = 12(2) in 3D at 
T = 1 and x = 5.3(3) in 4D at T = 1/ 0.6 which are both much larger than the mean-field 
value of X = 3. 

By looking for reasonable alternatives we realized that for many other systems the 
statistics of extremes as pioneered by Frechet, Fisher and Tippert, and von Mises"'^ '^^ has 
let to a good ansatz with universal properties^®'^^. It is based on a standard result^**'^^, due 
to Fisher and Tippert, Kawata, and Smirnow, for the universal distribution of the first, sec- 
ond, third, . . . smallest of a set of N independent identically distributed random numbers. 
For an appropriate, exponential decay of the random number distribution their probability 
densities are given by the Gumbel form 

Ux) = Ca exp[a (a;-e^)] , (21) 

in the limit of large A^. The exponent a takes the values a, = 1,2.3,..., corresponding, 
respectively, to the first, second, third, . . . smallest random number of the set, x is a scaling 
variable which shifts the maximum value of the probabiUty density to zero, and Ca is 
a normaUzation constant. For certain spin-glass systems the possible relevance of this 
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Figure 8. Re.scaled overlap probability densities for the EAI spin-glass model on L'^ lattices at the critical tem- 
perature on a linear (left) and expanded logarithmic (right) scale. In the lower part of the plots the deviation 
Pjg(<;') — Pgj(g') ± APlg{q') of some L = 16 data from the moditied Gumbel tit (22) is shown (with an 
unimportant offset added in order to fit inside the figure). 



universal distribution has been pointed out by Bouchaud and Mezard^^, and for instance 
also in applications to the 2D XY model in the spin-wave approximation^^'^^ the Gumbel 
ansatz (21) fitted well, albeit with a modified value of a = 7r/2. 

In our case we set x = b{q' — g^ax) modified the first x on the r.h.s. of (21) to 
c tanh(x /c) , where c > is a constant, in order to reproduce the flattening of the densities 
towards q' = 0. Notice that in the tails of the densities, i.e. for large, positive x, this term is 
anyway subleading. A symmetric expression for P'{q') reflecting its q' —q' invariance 
is obtained by multiplying the above construction with its reflection about the g'' = axis, 

P'{q') = Cexp |a ctanh (^^ {q' - q'^,^)^ - e+^^'-w) | x (g' ^ -q') . (22) 

Of course, the important large-a; behaviour of Eq. (21) is not at aU affected by our manip- 
ulations. 

By fitting this ansatz to our 3D data we obtained final estimates of a = 0.448(40) for 
T = 1.14 and a = 0.446(37) for T = 1, respectively. For T = 1.14 our best fit is already 
included on the l.h.s. of Fig. 8. We see a good consistency between the data and the fit 
over a remarkably wide range of q' . Even more impressive is the excellent agreement in 
the tails of the densities. Taking the T = 1.14, L = 16 result at face value, we find^^ a 
very good fit over a remarkable range of 200/ In(lO) « 87 orders of magnitude! 



4 Summary and Conclusions 

Employing non-Boltzmann sampling with the multi-overlap MC algorithm we have inves- 
tigated the probability densities Pj{q) of the Parisi order parameter q. The free-energy 
barriers F^^ as defined in Eq. (13) turn out to be non-self-averaging. The logarithmic scal- 
ing ansatz (16) for the barriers at fixed values F of their cumulative distribution function 
is found to be favoured over the mean-field ansatz (15). Further, relevant barriers are still 
reflected in the autocorrelations of the multi-overlap algorithm. 
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The averaged densities exhibit a pronounced FSS collapse onto a common master curve 
even in the spin-glass phase. For the scaling of their tails towards = ±1 we find quaU- 
tative agreement with the decay law predicted by mean-field theory, but with an exponent 
X that is, in particular in 3D, much larger than theoretically expected. A much better fit 
over more than 80 orders of magnitude is obtained in 3D by using a modified Gumbel 
ansatz, rooted in extreme order statistics^'* "^'^. The detailed relationship between the EAI 
spin-glass model and extreme order statistics remains to be investigated, and it is certainly 
a challenge to extend the work of Bouchaud and Mezard^^ to the more involved scenarios 
of the replica theory. 
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